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We report the ab initio study of rare-earth adatoms (Gd) on an insulating surface. This surface is 
of interest because of previous studies by scanning tunnehng microscopy showing spin excitations of 
transition metal adatoms [1]. The present work is the first study of rare-earth spin-coupled adatoms, 
as well as the geometry effect of spin coupling, and the underlying mechanism of ferromagnetic cou- 
pling. The exchange coupling between Gd atoms on the surface is calculated to be antiferromagnetic 
in a linear geometry and ferromagnetic in a diagonal geometry, by considering their collinear spins 
and using the PBE+U exchange correlation. We also find the Gd dimers in these two geometries 
are similar to the nearest-neighbor (NN) and the next-NN Gd atoms in GdN bulk. We analyze how 
much direct exchange, superexchange, and RKKY interactons contribute to the exchange coupling 
for both geometries by additional first-principles calculations of related model systems. 



PACS numbers: 

Understanding the spin coupling at the nanoscale is 
important in scaling down magnetic devices, such as spin- 
tronics and quantum computing devices [2|. During the 
past years, it has been demonstrated that the scanning 
tunnelling microscope (STM) is a powerful tool to build 
and control individual nano magnetic structures, and has 
a great potential to being applied to the construction of 
nanoscale magnetic devices. Previous studies of engineer- 
ing individual magnetic atoms on surfaces include the 
antiferromagnetism of Mn chains (1 ~ 10 atoms) [l|, the 
anisotropy of a single Fe atom [3] and a Fe-Cu dimer [^ , 
Kondo effects under magnetic environments [5, 6], and 
a bistable atomic-scale Fe antiferromagnet that demon- 
strates the feasibility of dense nonvolatile storage of infor- 
mation at low temperature [7] , where all above magnetic 
systems are placed on top of a copper nitride island that 
serves as an insulating monolayer on a Cu(lOO) substrate. 

Before physicists started using the STM to manipu- 
late and couple magnetic atoms together, chemists have 
decades of history in synthesizing numerous species of 
molecules that carry giant spins, well-known as molecu- 
lar magnets [8]. The STM-engineered spins were found 
to form a surface molecular network with great similarity 
to molecular magnets. While the major attention to the 
field of molecular magnets is mainly focused on those con- 
sisting of transition- metal magnetic atoms, very few stud- 
ies are devoted to the rare-earth-based molecular mag- 
nets. Moreover, only magnetic anisotropy is studied for 
the rare-earth-based molecular magnets [9], but not the 
interatomic spin coupling within such molecules. A sim- 
ilar situation also happens to the STM-engineered spins; 
experimentalists have not tried to place rare earth atoms 
on the CuN/Cu(100) surfaces to see how their spins cou- 
ple to each other. 

Rare earth atoms have their magnetism primarily con- 
tributed from the / orbit als, and may behave quite differ- 
ently from those well-studied transition atoms when be- 



ing placed on the surface. Among the lanthanoid series, a 
free Gd atom has a half-filled 4/ shell such that it carries 
an atomic spin much larger than the transition atoms, 
and at the same time has an / shell of L = that is ex- 
pected to exhibit quite small magnetic anisotropy on sur- 
faces due to small spin-orbit interaction. The spin excita- 
tions of spin coupling and magnetic anisotropy can both 
contribute to the STM inelastic tunnelling spectroscopy. 
If both types of excitations exist at the same energy scale, 
it makes the inelastic tunnelling spectra difficult to be 
analyzed. On the other hand, with low anisotropy one 
would expect that Gd atoms on the CuN surface yield 
clean inelastic tunnelling spectra mainly from the inter- 
atomic coupling of their spins. Such an advantage would 
benefit future experimental studies following this first- 
principles investigation. 

With the goal of making nano magnets by atom ma- 
nipulation, one needs to couple atomic spins ferromag- 
netically on such a surface. Previous first-principles cal- 
culations show the antiferromagnetism originated from 
superexchange interactions along a linear path [10| . The 
Gd atoms also serve as a potential candidate for ferro- 
magnetic atomic-spin dimers, as we describe below. 

In this work, we perform first-principles calculations 
of Gd adatoms on the CuN surface. In some ways, the 
Gd atoms are similar to the previously-studied Mn atoms 
[3|, llO| when being deposited on the Cu sites of the CuN 
surface, i.e., the Gd's nearby N atoms break bounds with 
their neighboring Cu and form a "quasi" molecular struc- 
ture from the surface. However, the local structures of 
the Gd atoms on the CuN surface have a well-studied 
reference system, the GdN bulk. We build two different 
geometries of the Gd dimers on the surface: one has Gd 
atoms along the same N row, and the other along two 
perpendicular N rows. The two geometries mimic the 
coupling paths of the NN and next-NN Gd atoms of the 
GdN bulk, where the two paths in bulk have ferromag- 




Figure 1: The unit cell of a GdN bulk. The Gd-to-Gd arrows 
indicate the NN (diagonal, purple) and next-NN (linear, red) 
couplings. 
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Figure 2: The upper figure shows the top view of the relaxed 
structure of a linear Gd dimer on the CuN surface. The lower 
left (right) figure shows the top view of the initial (relaxed) 
structure of a diagonal Gd dimer on the CuN surface. Unit 
cells are marked by red rectangles. 



netic and antiferromagnetic couplings, respectively. We 
calculate the exchange couplings J of two arrangements 
of Gd2/CuN using first-principles PBE+U, and expect 
that one of the two types of surface Gd dimers will ex- 
hibit ferromagnetism and the other antiferromagnetism. 
In the STM experiments, a copper-nitride monolayer is 
built between a magnetic atom and the Cu(lOO) surface 
[l|, y-lsl to keep the atomic spin away from the screening 
of its underlying conduction electrons while permitting a 
sufficient amount of STM tunneling current for probing 
the spin excitations. To understand the magnetic prop- 
erties of Gd atoms on the CuN surface, we simulate a 
single Gd on this surface by first constructing a supercell 



of 5-layer Cu slabs plus 8 vacuum layers with the nitro- 
gen atoms snugging in-between the half of vacant sites 
and then placing each Gd atom atop the Cu site of the 
CuN surface within a 3 unitcell. We perform density- 
functional calculations in the all-electron full-potential 
linearized augmented plane wave (FLAPW) basis |11|. 
A naive local density approximation (LDA) or gener- 
alized gradient approximation (GGA), when being ap- 
plied to materials composed of rare-earth atoms, gener- 
ally yields / levels inconsistent with photoemission ex- 
periments, and needs to be fixed by adding extra on- 
site Coulomb repulsion to the exchange-correlation func- 
tional, the so-called DFT+U method. To determine the 
on-site Coulomb Uf and exchange Jf values of the Gd 4/ 
orbit als on the CuN surface, we revisit the GdN bulk sys- 
tem, which mimics very well the local structure of Gd on 
the CuN surface. The GdN bulk, being a ferromagnetic 
semiconductor, has been studied experimentally by pho- 
toemission and computationally by LSDA+U for its po- 
tential application in spintronics. Following the same way 
as the previous LSDA+U studies of GdN bulk [H, [III, 
we find that Uf = 6.7 eV and Jf = 0.7 eV yield the 
energy difference between the majority-spin Gd 4/ and 
N 2p states in best agreement with photoemission mea- 
surements 15|, ll6| . This set of Uf and Jf are used in our 
succeeding calculations of both the GdN bulk and the Gd 
dimer on the CuN surface. 

The spin couplings along the diagonal and the linear 
Gd-N-Gd paths of the GdN bulk are well-studied in the 
literature [13|, [l^, where the two Gd atoms along the 



paths are NN and next NN to each other, respectively, 
as shown in Fig. [H The two coupling paths between Gd 
atoms in a GdN bulk strongly suggest that there are two 
possible geometries of Gd dimers on the CuN surface: 
one has Gd-N-Gd along the N row, and the other in a 
right angle (two Gd atoms along the diagonal) . Previous 
GdN-bulk studies have concluded that the spin couplings 
between Gd atoms of a GdN bulk along the diagonal and 
the linear paths are ferromagnetic and antiferromagnetic, 
respectively. Therefore we expect the surface Gd dimers 
in two geometries to have spin couplings the same as 
their counterparts of similar geometries in the GdN bulk, 
i.e. diagonal (linear) being ferromagnetic (antiferromag- 
netic). 

We therefore arrange Gd atoms in those two geometries 
on the CuN surface, and optimize the crystal structures 
until the maximum force among all the atoms reduces 
to ^10 mRy/ao and 5 mRy/ao- The relaxed structures 
are shown in Fig. [21 and the dimer local geometries are 
quantitatively presented in Table HI It is interesting to 
notice that the diagonal Gd dimer relaxes its bond angle 
from 90° to 112°. This can be understood given that the 
diagonal Gd-to-Gd distance in GdN bulk is 3.52A, and 
the initial Gd-Gd distance on the surface is 2. 56 A, much 
shorter than 3.52A, so a relaxed Gd-Gd distance on the 
surface of 3. 64 A is rather reasonable. To determine the 




-4 -2 

Energy (eV) 



Figure 3: Partial density of states of a single Gd on the CuN 
surface. 



Gd spin on the CuN surface, we plot the calculated par- 
tial density of states (PDOS) of a single Gd on the CuN 
surface in Fig. [H One clearly sees that the 4/ majority 
spin states are all occupied and the minority states are 
all unoccupied, which implies a 4/^ configuration for Gd, 
a spin-7/2 configuration for its 4/ shell. In addition, the 
5d states are not occupied, and its rather small PDOS in 
the entire energy range indicates its delocalization out- 
side of the Gd atomic sphere, in contrast to a free Gd 
atom that carries a valence configuration 5(i^4/^6s^. By 
comparing the Gd2/CuN and the GdN bulk, we find that 
both systems have each Gd atom connecting to N atoms 
such that their Gd local structures are very similar to 
each other but are significantly different to that of a free 
atom. The local structure plays an important role of the 
spins of all three Gd-contained systems. 

To calculate the spin coupling J between Gd spins on 
the CuN surface, we take advantage of the correspon- 
dence between the collinear spins of a Heisenberg model 
and the magnetic moments of the real crystal surface of 
interest [10] . The Hamiltonian of a Heisenberg spin dimer 



H = JSi-S2. 



(1) 



The difference of energy expectation values SE between 
the parallel and antiparallel spins is related to the cou- 
pling J, for spin-^* atoms, as 



SE = JS^ - i-JS^) = 2JS^ 



(2) 



By calculating the total energies of the parallel- and 
antiparallel-spin configurations of a Gd dimer at the 
Cu site of a CuN surface, we obtain from (J2j) the ex- 
change coupling J to be 1.24meV for the linear dimer 
and — 1.25meV for the diagonal dimer. Our calculations 
obtain an antiferromagnetic coupling for the linear Gd 





Gd to Gd (A) 


Gd-N (A) 


ZGd-N-Gd 


J (meV) 


NN Gd 

in GdN bulk 


3.52 


2.49 


90° 


-1.51 


next NN Gd 
in GdN bulk 


4.98 


2.49 


180° 


1.09 


diagonal Gd dimer 
on CuN surface 


3.64 


2.19 


112° 


-1.25 


linear Gd dimer 
on CuN surface 


4.15 


2.24 


135° 


1.24 



Table I: Calculated Gd-to-Gd distances, Gd-N bond length, 
Gd-N-Gd bond angle, and the spin coupling J between Gd of 
four systems. 





J 


Jo. 


Jp 


Jd 


Js 


Jr 


diagonal 


-1.25 


-0.92 


-1.30 


-0.97 


-0.33 


0.05 


linear 


1.24 


-0.88 


0.66 


-1.46 


2.12 


0.58 



Table II: Calculated spin couplings in meV of the original Gd 
dimers on the CuN/Cu(100) surface J, the Gd dimers on a 
single CuN sheet J/3, and Ja the spin coupling for the case 
of the N in between the Gd substituted by a Ne atom. Also 
listed are direct exchange J^, superexchange Js, and RKKY 
Jr extracted from J, J^, and Jp. 



dimer and ferromagnetic for the diagonal dimer. One 
can also see in Table |T] a trend of the coupling J vary- 
ing with the angle ZGd-N-Gd. The two Gd atoms start 
with a ferromagnetic coupling at an exact right angle, be- 
comes slightly less ferromagnetic at ZGd-N-Gd = 112°, 
changes to antiferromagnetic at ZGd-N-Gd = 135°, and 
finally stays antiferromagnetic when along a straight line. 

We now turn our attention to the underlying mecha- 
nism of the ferromagnetism (antiferromagnetism) of the 
diagonal (linear) Gd dimer. There are three possible 
magnetic interactions that may couple the two Gd spins: 
direct exchange, superexchange, and the Ruder man- 
Kittel-Kasuya-Yosida (RKKY) interaction. To extract 
the three components out of the resultant coupling, we 
perform calculations of two alternative model systems. 
One is the original Gd dimer on the surfce with the in- 
between N atom replaced by a Ne atom, which effectively 
turns off the superexchange. The other is a Gd dimer on 
top of a single CuN monolayer, i.e., removing the underly- 
ing metallic Cu(lOO) slab, with the dimer and CuN sheet 
remaining, which has basically no RKKY. The atomic 
positions of these two alternative model systems exactly 
follow the original realistic surface, i.e. their structures 
are not relaxed, so that the magnetic coupling is the only 
difference among these systems. We then decompose the 
spin couplings of the three systems: the original Gd dimer 
on the CuN/Cu(100) surface J, the one with a Ne in be- 
tween two Gd atoms Joc^ and the Gd dimer on a single 
CuN sheet J^, into the contributions of direct exchange 





Figure 4: Calculated spin-density isosurfaces of diagonal (up- 
per) and linear (lower) Gd dimers on CuN in their ferromag- 
netic (left) and antiferromagnetic (right) configurations. Red 
stands for positive magnetization, and blue for negative. The 
white arrow in the stick-ball inset of each plot indicates the 
observation direction. The green and yellow solid circles de- 
note the positions of N and Cu, respectively. All spin densities 
are plotted at the magnitude of O.OOle/A^. A "*" symbol is 
labeled in the less stable spin configuration of each Gd dimer 
geometry. 



Jd^ superexchange J^, and RKKY J^. 
their relations as 

J ^ Jd -p Jg \ ^ri 
J OL ^ J d ~r Jr-) 

Jl3 = Jd-\- Js- 



We write down 



(3) 



It is a simple matter to solve for Jd, Js-, and Jr- The 
calculated spin couplings are listed in Table [III 

For the diagonal dimer, the insignificant differences 
among J, J ex, and Jg imply that the N atom in between 
and the underlying conduction electrons play minor roles 
in the spin coupling of the Gd dimers, while the direct 
wavefunction overlapping between the two Gd atoms ac- 
tually dominates. In fact, the obtained Jd, Js-, and Jr 
values of the diagonal dimer reflect the statement above, 
where the superexchange and RKKY are 34% and 5% 
of the direct exchange, respectively. Notably, we find a 
ferromagnetic superexchange. The linear dimer, on the 
other hand, has quite different components of the three 
types of magnetic couplings. Antiferromagnetic superex- 
change dominates, larger than the direct exchange, and 
there is also significant antiferromagnetic RKKY, at 27% 
of superexchange. When comparing the same types of 
magnetic interaction between the two dimer geometries, 
we find that the direct exchange of both dimers have the 
same sign and order of magnitude while their superex- 
change Js have opposite signs. Consequently, the mag- 
nitude of Js determines the sign of the dimers' total J. 



Figure 5: Calculated spin-density isosurfaces of the modeled 
Gd dimers on CuN with the in-between N replaced by Ne. 
Such dimers will be called Gd-Ne-Gd dimer later, and both 
are in their ferromagnetic configuration here, (a) and (b) 
The side and top views of the diagonal, (c) and (d) The side 
and top views of the linear. For all four plots, red stands 
for positive magnetization, and blue for negative. The white 
arrow in the stick-ball inset of each side-view plot indicates the 
observation direction. The green, yellow, and light blue solid 
circles denote the positions of N, Cu, and Ne, respectively. All 
spin densities are plotted at the magnitude of 0.0005e/A^. 



With the knowledge that the ferromagnetic coupling for 
the diagonal case is partly due to ferromagnetic superex- 
change, we revisit Table HI and interpret the trend of the 
coupling J varying with the angle ZGd-N-Gd, as a rough 
measure of the angle dependence of superexchange in 
these systems. The prediction of ferromagnetic superex- 
change at 90 degrees, varying through a sign change and 
ending up at antiferromagnetic for 180 degrees is well 
known in the chemistry community [18], and is due to 
different orbitals being involved with the superexchange 
hopping to the N at different angles. Pauli exclusion 
favors antiferromagnetic coupling only when the three 
atoms are more or less in a line. At sharper angles, the 
relative symmetries of the orbitals involved in the hop- 
ping through the N tend to favor ferromagnetic alignment 
of the spins. 

In order to obtain more physical insights for the 
Gd dimers, we plot their spin densities for both the 
parallel- and antiparallel-spin configurations in Fig. HJ 
The parallel-spin configuration of a diagonal dimer has 
its spin density forming one isosurface lobe, while the 
antiparallel-spin has two disconnected Gd lobes of oppo- 
site spins. In the ferromagnetic configuration, the inter- 
mediate N spin is partially enveloped by the Gd-dimer 
spin lobe, and is carrying an oppositely polarized spin. 
In the antiferromagnetic configuration, the N atom in be- 
tween becomes a magnetic dipole antisymmetrically po- 
larized by the two opposite Gd spins. We compare this 



spin density with that of a hnear Gd dimer on the same 
surface in Fig. HI We see that while the ferromagnetic di- 
agonal dimer has its two spin lobes connected with each 
other, the ferromagnetic linear dimer forms two disjoint 
ones. The linear dimer has an antiferromagnetic ground 
state, and the corresponding spin density has a nodal 
plane exactly in the middle of the two Gd. Direct spin 
interchange in the linear case is expected to be less strong 
because of the intervening Nitrogen, and compared to the 
"spin bonding" between the two Gd atoms in the diago- 
nal case, the latter implying a strong overlapping of their 
spin unpaired orbit als. However, from Table [III we see 
that the linear case has direct exchange 1.5 times that of 
the diagonal. 

In order to understand the physics underlying the di- 
rect exchange in these two configurations, we set out to 
isolate the direct exchange. The strong superexchange of 
the linear dimer results in an antiferromagnetic ground 
state, and the dimer has a fundamental nodal plane in the 
middle, which obscures the underlying direct exchange in 
the spin density. To eliminate the superexchange from 
the system, we study instead our model systems with N 
replaced by Ne. Both superexchange- free configurations 
exhibit a ferromagnetic ground state, which we attribute 
to the underlying direct exchange. We plot their ferro- 
magnetic spin densities in Fig. [5l The diagonal Gd-Ne- 
Gd dimer has a spin-density lobe very similar to that of 
a realistic Gd dimer, which is consistent with the domi- 
nant direct exchange seen in Table |TT1 The linear Gd-Ne- 
Gd dimer, on the other hand, provides new information 
from its shape. We see a negative spin lobe exactly in the 
middle, together with prolongation of positive polariza- 
tion perpendicular to the dimer, as viewed from the top. 
In both configurations we see the strong effect of orbital 
symmetry in the interactions. In the diagonal case there 
are no nodes, and the linear case we see there are two 
nodes, with a small maximum in the middle. This strong 
indication of orbital interaction is further confirmation 
that we are seeing direct exchange. 

In order to understand further the details of the in- 
terstitial spin density between the Gd atoms of the di- 
agonal dimer, we plot the spin-dependent partial den- 
sity of states (S-PDOS) in the entire interstitial region 
of the unit cell (Fig. [6^) where S-PDOS is defined as the 
minority-spin DOS subtracted from the majority-spin in 
this region. There is an isolated positive peak at £^ = — 9 
eV, which is exactly the energy of the majority-spin 4/ 
level. The highly oscillating S-PDOS within the energy 
range from £^ = — 8 to — 1 eV results in the spins of 
the majority- and minority-spin orbitals approximately 
cancelling. In contrast, from E = —1 eV to the Fermi 
energy, the S-PDOS shows a peak that dominates, la- 
belled "1^1 ". Di actually represents three equivalent near- 
energy peaks. There are also two other peaks in region 
A. In Fig. [6t, we plot the majority-spin probability den- 
sity of Di. We can see that the highest probability of 



CO 

o 

Q 
Q. 




-6 -4 

Energy (eV) 



40- 



w, 10- 
V) 

O 

Q 
Q. 



-10- 



-20- 



-30- 



4f^— Mi 




-10 



-6 -4 

Energy (eV) 







Figure 6: (a) Spin-dependent partial density of states (S- 
PDOS) in the entire interstitial region of a ferromagnetic di- 
agonal Gd-Ne-Gd dimer on CuN, defined as the majority-spin 
DOS minus the minority-spin in this spatial region. The plot- 
ted energy range is divided into three regions. Region A: the 
S-PDOS is mostly positive. The peak that dominates the S- 
PDOS is labelled as ''Di'\ Region B: the S-PDOS is highly os- 
cillating. Region C: a single positive peak at the energy of the 
majority-spin 4/ level. The inset shows the Gd 5d S-PDOS 
in the energy range A. (b) Same plots as (a) for the linear 
Gd-Ne-Gd dimer. The dominating peak is labelled as "Li". 
(c) The majority-spin probability densities of the orbitals Di 
(left) and Li (right), viewed along the same orientation as 
the side view of the spin density of the diagonal and linear 
Gd-Ne-Gd dimers in Fig. [S] respectively. The isosurfaces are 
plotted at the value of 0.02A~^. 



Di is mainly concentrated in the interstitial region in be- 
tween two Gd, and provides the major contribution to 
the interstitial spin density. Although we do not plot the 
other two near-energy peaks here, we find that they look 
very similar. In contrast, all other peaks, positive and 
negative, in region A have negligible weight between the 
two Gd, and represent the interstitial electrons in the 
bulk Cu. The interstitial contribution of Di is decom- 
posed into plane waves in the FLAPW basis. To trace 
the atomic configuration of the above interstitial spin be- 
tween the Gd, we calculate the projections of the Di or- 
bital to s, p and d symmetries within a Gd muffin-tin 
sphere. (There is no appreciable / DOS in the inter- 
stitial region.) The calculated muffin-tin projections to 
6s, 6p, and 5d are 17%, 13.5%, and 68.5% of Gd, re- 
spectively. We also plot the Gd 5d S-PDOS around the 
Di energy, and find that there is also a peak exactly at 
Di. Consequently we believe that the Gd 5d electrons 
are delocalized across the diagonal to enable the direct 
exchange. 

We also perform S-PDOS analysis for the linear Gd- 
Ne-Gd dimer, and plot it in Fig. ^jp. The linear-dimer 
S-PDOS can be analyzed in a manner similar to that of 
the diagonal. There are a number of peaks in region A. 
However, the peak that dominates the S-PDOS surround- 
ing the Gd is the one which we label "Li" in the figure. 
Unlike the diagonal Gd-Ne-Gd dimer, the Li peak has 
a wide width and does not have energy neighbors. The 
majority-spin probability is also different from Di, and 
is elongated on the surface perpendicularly to the Gd-Gd 
direction, very similarly to Fig. [SJi. As seen in Fig. [6h, 
we notice that the two Gd orbital lobes of Li, unlike the 
1^1, do not connect with each other, and instead have a 
node in between and an antibonding appearance. 

As calculated for Di, the muffin-tin projections of Li 
to Gd 6s, 6j9, and 5d are 33%, 3%, and 64% of Gd, re- 
spectively. The 5d character of Li within the Gd muffin- 
tin sphere is consistent with the bd S-PDOS plotted in 
the inset of Fig. [BJd. Therefore, we also believe that the 
Gd bd electrons, as in the diagonal case, are delocalized 
to enable the direct exchange of the linear dimer. The 
spreading of Li out of the Gd muffin- tin spheres into the 
interstitial region represents 54% of the total probability 
distribution, while Di plus its two near-energy neighbors 
represent only 43%. We see the linear configuration has 
a contribution roughly 1.3 times that of the diagonal, a 
ratio close to the 1.5 ratio between the direct exchange 
couplings of the two configurations shown in Table [III In 
fact, if we quantify the interstitial magnetic moments of 
the diagonal and linear Gd-Ne-Gd dimers, they are 1.53 
and 2.34, respectively, with a ratio 1.54, even closer to 
1.5. The absence and presence of a middle nodal plane 
in between the two geometries of the Gd-Ne-Gd dimers 
are very likely due to the symmetry of their (hybridized) 
delocalized Gd orbitals, a coupling we refer to generically 
as spin bonding. 
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Figure 7: (a) Calculated spin-density isosurfaces of a GdN 
bulk. The right plot is viewed along the direction perpendic- 
ular to a cube face, and has the isosurfaces of only the nearest 
atoms plotted for better visualization. The plotted region rep- 
resents a GdN conventional unit cell extended by 30% of the 
lattice constant, so that the spin isosurfaces of atoms at face 
centers and corners can be entirely plotted. Red stands for 
positive spin polarization, and blue for negative. Gd and N 
atoms have nearly spherical positive and negative isosurfaces, 
respectively, (b) Interstitial and Gd 5d PDOS of a ferro- 
magnetic GdN bulk. The main peaks of occupied majority- 
and minority-spin PDOS are labeled as Bup and Bdn, respec- 
tively, (c) Isosurface of the spin probability density of the 
PDOS peak at the value of 0.0002A-^ 



We now consider bulk GdN. A natural generalization of 
the spin bonding of the diagonal-dimer direct exchange 
is to ask whether the ferromagnetic NN coupling in a 
GdN bulk is, like the diagonal Gd dimer, also related to 
a spin bonding. For this we plot the spin density of the 
GdN bulk in Fig. [7^, where it can be clearly seen that 
there is no obvious spin density lobe connecting two Gd 
atoms, i.e., no spin bonding in a GdN bulk. In fact, 
a previous DFT study of GdN [12] already concludes 
that the ferromagnetism of the diagonal nearest neigh- 
bors most likely originates from the RKKY interaction, 
by analyzing the trends of exchange coupling of differ- 
ently strained GdN compared to all other gadolinium 
pnictides of larger-size anions. RKKY is quite prominent 
in Gd bulk compounds, compared to a pair of surface 
adatoms, because of the increased free-electron density 
of states contributed by all the Gd. As we will show in 
the next paragraph, the oscillating spin density in the 
interstitial region makes direct exchange much less likely 
in the bulk. 

To understand how Gd spins couple each other through 
their delocalized electrons in a GdN bulk, we plot both 
the interstitial and 5d PDOS. This is shown in Fig. [Tb. 
Within 5eV below the Fermi energy, the main PDOS 
peak around E = —1.2 eV has a significant interstitial 
component but only a minor 5d. We further plot both the 
majority- and minority-spin orbitals associated with this 
peak, shown in Fig. [7t. The isosurfaces of the orbitals 
spread all over the interstitial region. This can be inter- 
preted as the follows: The Gd 5d electrons in the GdN 
bulk join the conduction electron sea, and consequently 
they contribute to the RKKY interactions, but do not 
form a directional orbital interaction in between two Gd 
as in the surface dimer case. Although the ferromag- 
netism seems to be understood quite well by RKKY, the 
GdN-bulk study [12] does not exclude the possibility of 
a ferromagnetic contribution from the 90°-superexchange 
interaction which those authors note can be quite large in 
some oxides. Indeed, from our finding of superexchange 
which becomes ferromagnetic at an angle of 112 degrees, 
we predict that the 90-degree coupling in bulk material 
has at least some ferromagnetic superexchange as well. 

In summary, we have calculated the electronic struc- 
ture of coupled rare-earth (Gd) spins on a surface using 
the PBE+U exchange correlation, in the first study of 
rare-earth spin-coupled adatoms. The presence of Gd 
gives rise to rearrangement of the atomic structure that 
is quite different from that of a Mn atom [10]. The ge- 
ometry effect of the spin coupling is manifested by calcu- 
lating the exchange coupling between Gd atoms on the 
CuN surface and finding antiferromagnetic coupling in a 
linear geometry and ferromagnetic in a diagonal geome- 
try, showing that the sign of J can be tuned by different 
geometric arrangements. We also find the Gd dimer s 
in these two geometries have many similarities to the 
nearest-neighbor (NN) and the next-NN Gd atoms in a 



GdN bulk. 

The underlying physics of the dimer s' magnetism is 
studied by decomposing the magnetic couplings into the 
direct exchange, superexchange, and RKKY interaction, 
and the strength of the direct exchange is further pic- 
torially understood by the "spin bonding" between two 
Gd. The diagonal case has ferromagnetic contributions 
from both direct exchange and superexchange. The an- 
tiferromagnetism for the linear geometry is due to the 
predominance of superexchange for this configuration, 
notwithstanding a large ferromagnetic direct exchange. 
Superexchange is present in both geometries for the Gd 
dimer, ferromagnetic for the diagonal and antiferromag- 
netic for the linear. Even for the ferromagnetic coupling 
of the diagonal dimer, superexchange constitutes 27% of 
the total interaction. While the bulk GdN compound is 
basically dominated by the RKKY interactions, we find 
much smaller RKKY with the Gd dimers on a CuN sur- 
face. Our calculations also show that the Gd spin of 
these structures is 7/2, the same as that of GdN bulk, 
but different from a spin-4 free Gd atom with a valance 
configuration Af^bd^Gs'^. 
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